home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 23 / CU Amiga - Super CD-ROM 23 (June 1998).iso / CUCD / Graphics / STIMP_noise / source / ppmadapmedian / ppmadapmedian.c < prev    next >
Encoding:
C/C++ Source or Header  |  1998-01-31  |  8.8 KB  |  352 lines
  1.  
  2. /***********************************************************************/
  3. #define OP_NAME     "ppmadapmedian"
  4. #define VERSION     "1.05"
  5. #define DATE        "31.01.98" 
  6. #define AUTHOR      "Stefan Diener"
  7. /**********************************************************************/
  8.  
  9. #include <stdlib.h>
  10. #include <stdarg.h>
  11. #include <stdio.h>
  12. #include <string.h>
  13. #include <math.h>
  14. #include <sys/types.h>
  15.  
  16. #include <STIMP/ppm.c>
  17.  
  18. struct PPM_Info source, desti;
  19. static int m, n;
  20. static int edge=1, breite, mitte;
  21.  
  22. double mittel(unsigned char *quelle, int z, int s)
  23. {
  24.   int k, l, summe=0;
  25.  
  26.   /* Summe des Fensters */
  27.   for (k=-edge; k<=edge; k++)   /* Zeilen */
  28.     for (l=-edge; l<=edge; l++)   /* Spalten */
  29.       summe+=quelle[(z+k)*n+s+l];
  30.  
  31.   /* Summe/Anzahl berechnen */
  32.   return (double)summe/(double)(breite*breite);
  33. }
  34.  
  35. int median(unsigned char *quelle, int z, int s)
  36. {
  37.   static unsigned char bucket[256];
  38.   int k, l, summe=0;
  39.  
  40.   /* Buckets leeren */
  41.   for (k=0; k<256; k++) bucket[k]=0;
  42.  
  43.   /* Buckets fuellen */
  44.   for (k=-edge; k<=edge; k++)   /* Zeilen */
  45.     for (l=-edge; l<=edge; l++)   /* Spalten */
  46.       bucket[quelle[(z+k)*n+s+l]]++;
  47.  
  48.   /* Mitte suchen */
  49.   for (k=0; k<256; k++)
  50.   {
  51.     summe+=bucket[k];
  52.     if (summe>mitte) return k;
  53.   }
  54.  
  55.   /* durchgelaufen ??? */
  56.   return 255;
  57. }
  58.  
  59. void Do_It(void)
  60. {
  61.   char kandidat[27];
  62.   double resultat[27];
  63.   int i, j, k, l, a;
  64.   double r, g, b, h, s, u, v, rr, gg, bb, mini;
  65.   unsigned char matrix[3][3];
  66.   unsigned char *srcR, *srcG, *srcB, *dstR, *dstG, *dstB;
  67.  
  68.   srcR = source.redData;
  69.   srcG = source.greenData;
  70.   srcB = source.blueData;
  71.   dstR = desti.redData;
  72.   dstG = desti.greenData;
  73.   dstB = desti.blueData;
  74.   desti.maxval=source.maxval;
  75.  
  76.   for (i=edge; i<m-edge; i++)   /* Zeilen */
  77.   {
  78.     for (j=edge; j<n-edge; j++)   /* Spalten */
  79.     {
  80.       /* Mediane in jeder Komponente bestimmen */
  81.       matrix[0][0]=median(srcR, i, j);
  82.       matrix[1][1]=median(srcG, i, j);
  83.       matrix[2][2]=median(srcB, i, j);
  84.  
  85.       /* Median-Matrix fuellen */
  86.       for (k=-edge; k<=edge; k++)   /* Zeilen */
  87.         for (l=-edge; l<=edge; l++)   /* Spalten */
  88.         {
  89.           if (matrix[0][0]==srcR[(i+k)*n+j+l])
  90.           {
  91.             matrix[1][0]=srcG[(i+k)*n+j+l];
  92.             matrix[2][0]=srcB[(i+k)*n+j+l];
  93.           }
  94.           if (matrix[1][1]==srcG[(i+k)*n+j+l])
  95.           {
  96.             matrix[0][1]=srcR[(i+k)*n+j+l];
  97.             matrix[2][1]=srcB[(i+k)*n+j+l];
  98.           }
  99.           if (matrix[2][2]==srcB[(i+k)*n+j+l])
  100.           {
  101.             matrix[0][2]=srcR[(i+k)*n+j+l];
  102.             matrix[1][2]=srcG[(i+k)*n+j+l];
  103.           }
  104.         }
  105.  
  106.       /* Mittelwerte in jedem Kanal bilden */
  107.       r=mittel(srcR, i, j);
  108.       g=mittel(srcG, i, j);
  109.       b=mittel(srcB, i, j);
  110.  
  111.       /* RGB -> HSI: auf Singularitaeten pruefen */
  112.       if ((2.0*b-r-g==0.0) || ((r-g)*(r-g)+r*g+2.0*b*(b-r-g)<0.0))
  113.       {
  114.         /* Bedingungen nicht erfuellt => Originalpixel verwenden */
  115.         r=(double) srcR[i*n+j];
  116.         g=(double) srcG[i*n+j];
  117.         b=(double) srcB[i*n+j];
  118.  
  119.         /* Bedingungen nochmal pruefen */
  120.         if ((2.0*b-r-g==0.0) || ((r-g)*(r-g)+r*g+2.0*b*(b-r-g)<0.0))
  121.         {
  122.           /* auch nicht erfuellt => abbrechen */
  123.           r=0.0;  g=0.0;  b=0.0;
  124.         }
  125.       }
  126.  
  127.       if ((r!=0.0) || (g!=0.0) || (b!=0.0))
  128.       {
  129.         /* RGB -> HSI: Transformation */
  130.         u=2.0*b-r-g;
  131.         v=r-g;
  132.         h=atan2(v,u);
  133.         s=0.4082482905*sqrt(u*u+v*v);
  134.  
  135.         /* Listen leeren */
  136.         for (k=0; k<27; k++)
  137.         {
  138.           kandidat[k]=0;
  139.           resultat[k]=0.0;
  140.         }
  141.  
  142.         /* H-Transformation fuer alle 27 Kandidaten berechnen */
  143.         for (k=0; k<3; k++)   /* 1. Zeile */
  144.           for (l=0; l<3; l++)   /* 2. Zeile */
  145.             for (a=0; a<3; a++)   /* 3. Zeile */
  146.             {
  147.               rr=(double) matrix[0][k];
  148.               gg=(double) matrix[1][l];
  149.               bb=(double) matrix[2][a];
  150.  
  151.               /* Bedingungen ueberpruefen */
  152.               if ((2.0*bb-rr-gg!=0.0) && ((rr-gg)*(rr-gg)+rr*gg+2.0*bb*(bb-rr-gg)>=0.0))
  153.               {
  154.                 u=2.0*bb-rr-gg;
  155.                 v=rr-gg;
  156.                 resultat[a+l*3+k*9]=fabs(atan2(v,u)-h);
  157.               }
  158.               else resultat[a+l*3+k*9]=1000000.0;
  159.             }
  160.  
  161.         /* Minimum suchen */
  162.         mini=resultat[0];
  163.         for (k=1; k<27; k++)
  164.           if (resultat[k]<mini) mini=resultat[k];
  165.  
  166.         /* Kandidaten auswaehlen */
  167.         l=0;
  168.         for (k=0; k<27; k++)
  169.           if (resultat[k]==mini)
  170.           {
  171.             kandidat[k]=1;
  172.             l++;
  173.           }
  174.  
  175.         if (l>1)
  176.         {
  177.           /* S-Transformation fuer restliche Kandidaten berechnen */
  178.           for (k=0; k<3; k++)   /* 1. Zeile */
  179.             for (l=0; l<3; l++)   /* 2. Zeile */
  180.               for (a=0; a<3; a++)   /* 3. Zeile */
  181.                 if (kandidat[a+l*3+k*9]==1)
  182.                 {
  183.                   rr=(double) matrix[0][k];
  184.                   gg=(double) matrix[1][l];
  185.                   bb=(double) matrix[2][a];
  186.  
  187.                   u=2.0*bb-rr-gg;
  188.                   v=rr-gg;
  189.                   resultat[a+l*3+k*9]=fabs(0.4082482905*sqrt(u*u+v*v)-s);
  190.                 }
  191.  
  192.           /* Minimum suchen */
  193.           mini=1000000.0;
  194.           for (k=0; k<27; k++)
  195.             if ((kandidat[k]==1) && (resultat[k]<mini)) mini=resultat[k];
  196.  
  197.           /* Kandidaten auswaehlen */
  198.           l=0;
  199.           for (k=0; k<27; k++)
  200.             if ((kandidat[k]==1) && (resultat[k]==mini)) l++;
  201.            else kandidat[k]=0;
  202.  
  203.            if (l>1)
  204.            {
  205.             /* I-Transformation fuer restliche Kandidaten berechnen */
  206.             for (k=0; k<3; k++)   /* 1. Zeile */
  207.               for (l=0; l<3; l++)   /* 2. Zeile */
  208.                 for (a=0; a<3; a++)   /* 3. Zeile */
  209.                   if (kandidat[a+l*3+k*9]==1)
  210.                   {
  211.                     rr=(double) matrix[0][k];
  212.                     gg=(double) matrix[1][l];
  213.                     bb=(double) matrix[2][a];
  214.  
  215.                     resultat[a+l*3+k*9]=(rr+gg+bb)/3.0;
  216.                   }
  217.  
  218.             /* Maximum suchen */
  219.             mini=0.0;   /* Maximum heisst ausnahmsweise mini */
  220.             for (k=0; k<27; k++)
  221.               if ((kandidat[k]==1) && (resultat[k]>mini)) mini=resultat[k];
  222.  
  223.             /* Kandidaten auswaehlen */
  224.             for (k=0; k<27; k++)
  225.               if ((kandidat[k]==0) || (resultat[k]!=mini)) kandidat[k]=0;
  226.            }
  227.         }
  228.  
  229.         /* es bleibt nur ein Kandidat uebrig */
  230.         /* den Kandidaten suchen */
  231.         for (k=0; k<27; k++)
  232.           if (kandidat[k]==1)
  233.           {
  234.             l=k;
  235.             k=27;
  236.           }
  237.  
  238.         /* Daten aus Median-Matrix ins Zielbild kopieren */
  239.         k=l/9;  l=l%9;
  240.         dstR[(i-edge)*(n-2*edge)+j-edge]=matrix[0][k];
  241.         k=l/3;  l=l%3;
  242.         dstG[(i-edge)*(n-2*edge)+j-edge]=matrix[1][k];
  243.         dstB[(i-edge)*(n-2*edge)+j-edge]=matrix[2][l];
  244.       }
  245.       else   /* Singularitaet */
  246.       {
  247.         /* Quellpixel einfach in Zielbild kopieren */
  248.         dstR[(i-edge)*(n-2*edge)+j-edge]=srcR[i*n+j];
  249.         dstG[(i-edge)*(n-2*edge)+j-edge]=srcG[i*n+j];
  250.         dstB[(i-edge)*(n-2*edge)+j-edge]=srcB[i*n+j];
  251.       }
  252.     }
  253.   }
  254. }
  255.  
  256. int main(int argc, char **argv)
  257. /* Hauptprogramm */
  258. {
  259.   int i;
  260.  
  261.   /* offizielle Begruessung */
  262.   PrintOpening(argc,argv);
  263.  
  264.   /* Uebergebene Parameter auswerten */
  265.   for (i=1; i<argc; i++)
  266.   {
  267.     if ((argv[i][0]=='-') && argv[i][1])
  268.     {
  269.       switch (argv[i][1])
  270.       {
  271.         case '3': edge=1;
  272.                       break;
  273.  
  274.         case '5': edge=2;
  275.                       break;
  276.  
  277.         case '7': edge=3;
  278.                       break;
  279.  
  280.         case '9': edge=4;
  281.                       break;
  282.  
  283.         case 'v': beVerbose=FALSE;
  284.                       break;
  285.  
  286.         default: PrintMessage("Unknown parameter: %s", argv[i]);
  287.                      Hilfe();
  288.                      exit(-1);
  289.                      break;
  290.       }
  291.     }
  292.  
  293.     if (argv[i][0]=='+')
  294.     {
  295.       switch (argv[i][1])
  296.       {
  297.         case 'v': beVerbose=TRUE;
  298.                       break;
  299.  
  300.         default: PrintMessage("Unknown parameter: %s", argv[i]);
  301.                      Hilfe();
  302.                      exit(-1);
  303.                      break;
  304.       }
  305.     }
  306.   }
  307.  
  308.   /* Mindestzahl der Argumente pruefen */
  309.   if (argc<3)
  310.   {
  311.     PrintMessage("Wrong number of arguments !");
  312.     Hilfe();
  313.     exit(-1);
  314.   }
  315.  
  316.   /* Anzahl der Dateinamen überprüfen */
  317.   if (FilenameCount(argc, argv)!=2)
  318.   {
  319.     PrintMessage("Wrong number of file names !");
  320.     Hilfe();
  321.     exit(-1);
  322.   }
  323.  
  324.   if (ReadPPMFile(GetFilename(1,argc,argv),&source)==0)
  325.   {
  326.     m=source.height;
  327.     n=source.width;
  328.  
  329.     if ((m<=2*edge) || (n<=2*edge)) PrintMessage("The source picture is too small !!!");
  330.     else
  331.     {
  332.       if (CreatePPMArray(m-2*edge,n-2*edge,&desti)==0)
  333.       {
  334.         PrintMessage("Working ...");
  335.  
  336.         breite=2*edge+1;
  337.         mitte=2*edge*(edge+1);
  338.         Do_It();
  339.  
  340.         WritePPMFile(GetFilename(2,argc,argv),&desti);
  341.         FreePPMArray(&desti);
  342.       }
  343.     }
  344.  
  345.     FreePPMArray(&source);
  346.   }
  347.  
  348.   PrintClosing();
  349.   exit(0);
  350. }
  351.  
  352.